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We investigate the zero-temperature excitation spectrum of two-dimensional soft-core bosons for 
a wide range of parameters and across the phase transition from a superfluid to a supersolid state. 
Based on mean field calculations and recent Quantum Monte Carlo results, we demonstrate the 
applicability of the Bogoliubov-de Gennes equations, even at high interaction strengths where the 
system forms an insulating cluster crystal. Interestingly, our study reveals that the maximum 
energy of the longitudinal phonon band in the supersolid phase connects to the maxon energy of 
the superfluid at the phase transition. 
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A supersolid is a phase of matter that simultaneously 
accommodates diagonal as well as off-diagonal long-range 
order, which means that particles self- assemble into a 
rigid, regular crystal but at the same time can flow super- 
fiuidly through the formed solid. More than forty years 
ago, this peculiar state has been conjectured to emerge 
in pressurized solid Helium pj[3] , which led to an intense 
search for super solidity in such systems [4 j. In 2004, ex- 
perimental evidence for superfluidity in toroidal oscillator 
measurements [5j [6] has greatly revived interest in super- 
solid Helium and initiated a recent surge of activity on 
this problem. Yet, theoretical work has not reached a 
general consensus [7HTT] as to whether solid 4 He can dis- 
play superfluidity, and very recent experiments [12] are 
casting doubt on the original interpretation of the mea- 
surements [5j[6]. 

At the same time, ultracold atomic gases have emerged 
as a promising alternative platform to realize and study 
continuous-space supersolids in an unambiguous and con- 
trolled fashion. Recent work has demonstrated that long- 
range soft-core interactions between bosonic atoms can 
be engineered via optical coupling to highly-lying elec- 
tronic Rydberg states [T3H20] or via light-induced inter- 
actions in an optical cavity [21]. Such soft-core inter- 
actions can give rise to so-called cluster-solids [22] and 
cluster-supersolids [13] [15], i.e. crystalline arrangements 
of atomic clusters or droplets where superfluidity can 
arise from particle-hopping between the self-assembled 
droplets. 

On a more formal level, supersolidity can be under- 
stood in terms of the simultaneous breaking of funda- 
mentally different symmetries: (i) the breaking of trans- 
lational symmetry responsible for the crystalline ordering 
and {%%) the breaking a global gauge symmetry that en- 
ables long-range phase coherence and thereby superfluid- 
ity of the system. A direct consequence is the emergence 
of Goldstone bosons, i.e gapless modes in the excitation 
spectrum for each of the broken symmetries. Therefore, 
the excitation spectrum provides a powerful experimen- 
tal way to probe supersolidity, e.g. via Bragg scatter- 
ing [23j [24] and has recently attracted considerable the- 
oretical interest [25-31 . Most recently, Quantum Monte 
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FIG. 1: Superfluid fraction of two dimensional Bosons as 
a function of the dimensionless interaction strength a = 
mpVoRo/h 2 , that combines the particle mass m and den- 
sity p, as well as the strength Vb and range Ro of the 
soft-core interaction potential. The points show results of 
Monte Carlo simulations for 320 particles with a density 
of RqP = 4.4. The superfluid fraction drops discontinu- 
ously at a w 13.4 (dashed vertical line), marking a first- 
order superfluid-supersolid phase transition, in close agree- 
ment with the mean field prediction of a — 12.7 (solid verti- 
cal line). Around a ~ 22 quantum fluctuations destroy global 
phase coherence and superfluidity and the systems enters an 
insulating phase. The insets show PIMC snapshots for dif- 
ferent a, illustrating the particle density profile in the three 
different phases. 



Carlo [32] simulations have been used to determine the 
dynamical structure factor of two-dimensional soft-core 
Bosons [25], while mean field approaches have been ap- 
plied in numerous works [26l [27j I33H36] . aiming at a 
simplified description of the zero-temperature physics of 
these systems. 

Here we present a thorough comparison between these 
two approaches and show that the excitation spectrum 
can be accurately described by the Bogoliubov-de Gennes 
(BdG) equations. Surprisingly, we find quantitative 
agreement with Monte Carlo simulations, not only in the 
superfluid and supersolid phases but also in the strongly 



interacting regime, where superfluidity is destroyed by 
quantum fluctuations and mean field approaches are gen- 
erally expected to fail. This finding should prove valu- 
able for future work, as it justifies the use of more effi- 
cient mean field calculations, which moreover enable in- 
vestigations of dynamical processes, such as relaxation 
phenomena or externally driven systems. Here, we ex- 
ploit this fact to study the excitation spectrum for a wide 
range of parameters, and find unexpected features at the 
superfluid-supersolid phase transition. 
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I. GROUND STATE PROPERTIES 

We consider an ensemble of Bosons confined to two di- 
mensions with with mass m and positions q i5 as described 
by the Hamiltonian 



H 



E--v 



2m 



i<j 



(1) 



As a prototype example for soft-core interactions we 
chose a simple step function potential V(r) = VqO (Rq — 
g), where 0(r) denotes the Heavyside function and Vo 
and Rq define the strength and range of the interaction 
potential. Upon scaling lengths by Ro and energies by 
h 2 /mR^ the zero temperature physics, determined by 
eq.Q, depends only two the dimensionless parameters: 
an effective interaction strength a' = VomR^/% 2 and 
the dimensionless density pR^. At the mean field level, 
this set of parameters can be further reduced. In this 
case, the system dynamics is described by a non-local 
Gross-Pitaevskii equation (GPE), which, in terms of the 
rescaled variables, can be written as 



id t tp(r) 
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dr'C/(r-r')|^(r')| 2 



^(r) , (2) 



where r = q/i?o> U(r) = 0(1 — r), and a = mpV^R^jh 2 
is a dimensionless interaction strength that solely deter- 
mines the dynamics and ground state properties. Such a 
mean field treatment is expected to be valid in the limit of 
weak interactions, reached by decreasing a' and increas- 
ing the density p such that the product a' p stays finite. 
In this regime the zero temperature physics is determined 
only by the effective interaction strength a, which greatly 
simplifies the analysis of the underlying phase diagram. 

For small a the system is in a homogenous superfluid 
phase whose energy ira follows directly from eq.Q. We 
describe the modulated supersolid state by a variational 
wave function that is composed of localized Gaussians, 
arranged on a triangular lattice. Their dispersion a and 
the lattice constant a is obtained by minimizing the total 
energy. This simple analysis shows that density modula- 
tions become energetically favorable for a > 12.65, mark- 
ing a first order phase transition to a cluster supersolid 
state composed of small superfluid droplets p~3j [14] [32] . 
The droplets become more localized as the interaction 



FIG. 2: Excitation spectrum in the superfluid phase accord- 
ing to eq.Q (line) compared to the PIMC data (circles) of 
Ref. 25 . 



strength increases where both their size as well as the 
lattice spacing decreases from a = 0.39 and a = 1.51 at 
the phase transition to a = 0.22 and a = 1.40 for a = 40. 

In order to test these predictions we additionally per- 
formed path integral Monte Carlo (PIMC) simulations 
based on the continuous-space Worm algorithm [37] [38] . 
Fig{l] shows the obtained superfluid fraction as a function 
of a for a total number of 320 particles with a density 
of pRq = 4.4. One finds a first order phase transition at 
a ~ 13.4, signaled by a discontinuous drop of the super- 
fluid fraction from unity in the homogenous superfluid 
phase to ~ 0.6 at the supersolid transition. The location 
of the phase transition is in remarkable agreement with 
the above meanfield prediction, considering that there is 
only an average number of ~ 5.5 particles per droplet 
[43] . Also consistent with the above discussion, a fur- 
ther increase of a leads to a stronger localization of the 
droplets, and, hence, a drop of the superfluid fraction [3]. 
Around a « 22 the particle density between the droplets 
decreases to a point where quantum fluctuations destroy 
phase coherence between individual clusters such that the 
system enters an insulating crystal of superfluid droplets 
without long range off diagonal order [15 . While this 
transition can evidently not be captured by eq.Q, the 
meanfield theory nevertheless yields an accurate descrip- 
tion of the excitation spectrum in the insulating phase, 
as discussed below. 



II. EXCITATIONS 

The excitation spectrum is obtained by expanding the 
wavefunction, ip = ^o+c^, around the groundstate. Sub- 
stituting small perturbations of the form 



[u(r)e 



-iujt 



v*(r)e 



iut~\ 



(3) 



into the GPE Q yields to leading order in 5ij) the fa- 
miliar BdG equations for the Bogoliubov modes u and v. 
We solve these equations in real space by expanding the 
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FIG. 3: (color online) Mean field spectra (lines) at a = 16.93 (a) and a = 30.62 (b) obtained from eqs.|5} along the three 
symmetry directions of the Brillouin zone [see inset of panel (b)]. The symbols represent the PIMC data of Ref. [25] for 
longitudinal excitations computed along the direction r — M — V in the first two Brillouin zones. 



modes into Bloch waves 

u(r) = u k (r)e ikr 

v(r) = « k (r)e jk - r , (4) 

where the functions u\^(r) and v\^(r) obey the trans- 
lational symmetry of the underlying groundstate, i.e. 
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for and ^k, where A(r) = J dr f U(r f — r) | / 0(r , )| 2 . Re- 
cently, the excitation spectrum has been investigated by 
calculating the dynamical structure factor via PIMC sim- 
ulations [25 , using the so-called GIFT approach [39] . 
In the following, these first-principle results will be em- 
ployed to assess the predictive power of eq.([5|. 

In the homogeneous supernuid phase, eqs.Q can be 
solved analytically and yield the familiar Bogoliubov 
spectrum 

w = (y + 2atW), (6) 

where J7(k) = 2nJi(k)/k and J\ denotes the Bessel func- 
tion of the first kind. For a > 5.03 the spectrum develops 
a roton-maxon structure and roton softening occurs at 



a continuous translational symmetry in the supernuid 
phase and a discrete triangular lattice periodicity in the 
supersolid phase. This ansatz leads to the following set 
of equations 



- r>(rV k ' (r '~ r V(r') - «k(r')] = ™ k (r) 
-r')^(r')e ik -( r '- r )[u k (r')-t; k (r')] = wt; k (r) (5) 



a = 14.74, preceded by the supersolid phase transition 
at a = 12.7. Fig. [2] illustrates the roton-maxon spec- 
trum of the supernuid phase for an interaction strength 
of a = 11.86, i.e. slightly below the superfluid-supersolid 
phase transition. Eq.([6| is in excellent agreement with 
the numerical PIMC results of Ref. [25] , indicating that 
the system can indeed be described as a weakly coupled 
fluid. 

In the supersolid phase, a reliable calculation of the 
excitation spectrum requires accurate knowledge of the 
ground state and its chemical potential. To this end 
we iterate the time-independent GPE [40], fiipofjc) = 
[— V 2 /2+aA(r)]'0o(r), starting from our optimized varia- 
tional wave function and employing the same grid used to 
solve eqs.(|5|. In Fig. ^ we show the obtained spectrum 
of low-energy excitations for two values of the interac- 



^ V 2 + aA(r) - f?j u k (r) + ai/j(r) J drU(r 

• aip(r) J dr'U(r 



:V 2 +aA(r)-/i k(r) 



4 



15 



£10 

Oh 



bJO 
CD 

§20 

$=! 

.2 

1 10 







^^^^^ 





lb) 


1 1 1 1 1 












1 


i 






1 







10 20 30 

interaction strength, a 



40 



FIG. 4: (a) Linear dispersion of the gapless modes and (b) 
their energy at distinct momenta as a function of a. In the 
superfluid phase (a < 12.7) panel (b) shows the maxon energy 
(line) and the excitation energies (circles) of the three gap- 
less modes at the M-point (ujm) and of the longitudinal and 
transverse phonon band at the if -point (ojk) in the supersolid 
phase (see Fig|3]). 



tion strength that lie in the supersolid (a = 16.93) and 
in the insulating crystal phase (a = 30.62), respectively. 
The figure shows the excitation energies along the three 
symmetry axes of the Brillouin zone corresponding to 
the underlying triangular lattice. We find three gapless 
bands, i.e. three Goldstone modes reflecting the sym- 
metries that are broken in the supersolid phase [26]. In 
addition to the "superfluid band" due to the breaking of 
global gauge symmetry, there are two bands correspond- 
ing to longitudinal and transverse phonon excitations of 
the two-dimensional lattice. While the latter were not 
accessible by the PIMC calculations of Ref. [25] , we find 
good agreement for the two longitudinal modes in the su- 
persolid phase (a = 16.92). Somewhat surprisingly, even 
in the insulating phase, eq.([5| yields excellent agreement 
for the longitudinal phonon mode, despite its evident in- 
ability to describe the break-down of global superfluid- 
ity. This indicates that each individual droplet maintains 
a high condensate fraction despite the apparent lack of 
global phase coherence between the crystalline ordered 
droplets (see Figfl]). 

Having demonstrated the accuracy of the mean field 
description we can now exploit its efficiency to study 
the mode structure over a wider range of interac- 
tion strengths and, in particular, across the superfluid- 
supersolid phase transition. The results are summa- 
rized in Fig|4j which shows the linear dispersion v = 
duj/dk\k=o as well as the excitation frequency at several 
distinct momenta. As expected for a first order phase 
transition, the linear dispersion exhibits a discontinuous 
jump (FigJI^i), reflecting the sudden onset of finite den- 



sity modulations at the phase transition. On the other 
hand, the excitation energy shows a different and rather 
unexpected behavior. Figj4]3 shows the a-dependence of 
the maxon energy in the superfluid phase along with the 
excitation energy of the three gapless bands at the M- 
symmetry-point of the Brillouin zone (see Figj3|. At this 
symmetry point the longitudinal phonon mode provides 
the maximum energy in the entire reciprocal space. Inter- 
estingly, we observe that the maxon energy of the super- 
fluid merges into the latter at the superfluid-supersolid 
phase transition, i.e. does not show a discontinuous 
jump. Remarkably, at the phase transition, this energy 
also coincides with the excitation energy at the other 
symmetry point (if -point) where, in addition, the longi- 
tudinal and transverse phonon modes are degenerate. 



III. CONCLUSION 

Based on PIMC simulations and meanfield calcula- 
tions, we studied the zero-temperature physics of two- 
dimensional soft-core bosons and, in particular, their ex- 
citation spectrum at the phase transition to a supersolid 
droplet phase. We find remarkably good agreement be- 
tween the two approaches, even at densities for which 
the average particle number per droplet is as low as 
~ 5.5 and for interactions that are sufficiently strong to 
cause a breakdown of superfluidity. The close match be- 
tween both methods not only demonstrates the predic- 
tive power of the mean field approach but also attests to 
the accuracy of the GIFT method [39 for extracting dy- 
namical properties from Quantum Monte Carlo simula- 
tions. Our spectra are consistent with recent calculations 
of Refs. [26 , 27 , but differ qualitatively from the results 
of Ref. [41] , where only one gapless mode has been found 
in the supersolid phase of soft-core dipoles. 

A careful scan of the particle interaction revealed that 
the maxon energy of the superfluid phase merges into the 
energy of longitudinal phonons at both reciprocal sym- 
metry points of the supersolid phase. This degeneracy 
at the phase transition may suggest that the maxon part 
of the superfluid excitations plays a more significant role 
than thus far anticipated - an implication which calls 
for further investigation. Future studies of other types 
of soft-core interactions [13j [15j [21] will clarify whether 
this behavior is of generic nature or a mere consequence 
of the step-function potential considered in this work. 
Moreover, PIMC simulations [25 would allow to eluci- 
date the effects of correlations, i.e. to address the ques- 
tion whether the found connection between maxon and 
phonon excitations persists in the strong coupling regime 
beyond the validity of the BdG eqs.([5|. Along these lines, 
the behavior of Goldstone modes for related scenarios, 
e.g., at the crystallization point of dipolar systems [42] 
or at the superfluid-solid transition of Helium [11 , sug- 
gests itself as an interesting question for future studies. 
Finally, recent measurements that demonstrated super- 
solidity of ultracold atoms in an optical cavity via Bragg 
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scattering [21] would provide a viable experimental way 
to investigate this question for yet another important 
type of soft-core interactions. 
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